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Toward large eddy simulation 
of turbulent flow over an airfoil 

By Haecheon Choi 1 


1. Motivation and objectives 

The flow field over an airfoil contains several distinct flow characteristics, e.g. 
laminar, transitional, turbulent boundary layer flow, flow separation, unstable free- 
shear layers, and a wake. This diversity of flow regimes taxes the presently available 
Reynolds averaged turbulence models. Such models jure generally tuned to predict a 
particular flow regime, and adjustments are necessary for the prediction of a different 
flow regime. Similjir difficulties are likely to emerge when the large eddy simulation 
technique is applied with the widely used Smagorinsky model. This model has not 
been successful in correctly representing different turbulent flow fields with a single 
universal constant and has an incorrect near-wall behavior. 

Germano et al. (1991) and Ghosal, Lund & Moin (1992) have developed a new 
subgrid-scale model, the dynamic model, which is very promising in alleviating many 
of the persistent inadequacies of the Smagorinsky model: the model coefficient is 
computed dynamically as the calculation progresses rather than input a priori. 
The model has been remarkably successful in prediction of several turbulent and 
transitional flows. 

We plan to simulate turbulent flow over a “2D” airfoil using the large eddy sim- 
ulation technique. Our primary objective is to assess the performance of the newly 
developed dynamic subgrid-scale model for computation of complex flows about air- 
craft components and to compare the results with those obtained using the Reynolds 
average approach and experiments. The present computation represents the first 
application of large eddy simulation to a flow of aeronautical interest and a key 
demonstration of the capabilities of the large eddy simulation technique. 

2. Accomplishments 

2.1 Code modification 

For the simulation of the flow over an airfoil, we have modified the direct numer- 
ical simulation code of Choi, Moin it Kim (1993), which was successfully applied 
to simulation of turbulent flow over longitudinal riblets. The code uses generalized 
coordinates in two dimensions and a Cartesian coordinate in the third direction. 
The code is based on the unsteady, three-dimensional, incompressible Navier-Stokes 
equations along with the equation of continuity. The fully implicit time integra- 
tion scheme, which uses an approximate factorization method in conjunction with 
Newton-iterative scheme, allows rather large time-steps mid, therefore, reduces the 
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computational effort for solving the Poisson equation. The code makes use of a 
multi-grid method in order to accelerate the convergence rate of the Poisson solver. 

The convective outflow boundary condition is used at the flow exit, and the 
unsteady “turbulent” inflow condition is similar to that developed by Lee et al 
(1992). No-stress or uniform velocity boundary conditions are applied at far-field 
boundaries far from the airfoil surface. Periodic boundary condition is used in the 
spanwise direction. Arbitrary velocity boundary condition on the surface of the 
airfoil is also incorporated in the code for future flow control applications. 

The C-type mesh configuration is used in this study (Thompson et al. 1985). 
This grid is better than the O-type grid for flows with massive separation because 
the wake region after bluff bodies can be more adequately resolved with the C-type 
grid. The boundary condition along the branch-cut in the O-type grid is treated 
implicitly , i.e. periodic boundary conditions are applied along the branch cut, which 
makes the code significantly more complicated. The commonly used explicit inter- 
polation technique using adjacent grids near the branch cut deteriorates the time 
accuracy due to the time-lagging nature of the explicit method. However, the im- 
plicit boundary condition along the branch cut used in this study maintains the 
overall time accuracy. The use of explicit interpolation method along the branch 
cut does not accurately predict the Strouhal number for the flow over a cylinder, 
while the implicit boundary condition predicts the Strouhal number very well com- 
pared with the experimental results (see section 2.2). 

The development of the dynamic subgrid-scale model (Germano et al. 1991; 
Ghosal, Lund &: Moin 1992) is the key inducement for the undertaking of the 
present computation. The dynamic subgrid-scale model is incorporated in the code 
(in generalized coordinates). 

2.2 Code verification 

In order to verify the code, we first applied it to the laminar flow over a circular 
cylinder at the Reynolds number, Re = u^d/v = 100, where tioo is the upstream ve- 
locity, d is the cylinder diameter, and v is the kinematic viscosity. At this Reynolds 
number, the flow is periodic, and experimental data are available for comparison. 

The flow domain and grid system are shown in figure 1. The flow domain cov- 
ers about 25 and 20 cylinder diameters in the streamwise and normal directions, 
respectively. A uniform velocity field, u = iioo and v = 0, is prescribed at t = 0. 
Here, u and v are the velocity components in the streamwise (x) and normal (y) 
directions, respectively. An initial random disturbance with magnitude of 0.01 Uqo 
is imposed at the flow field in order to induce vortex shedding behind the cylinder. 

Figure 2 shows the time history of the streamwise velocity at several points behind 
the cylinder. After a transient period ( tu^jd ~ 0 ~ 100), the flow behind the 
cylinder shows a periodic behavior. The Strouhal number is calculated from figure 
2, 5 = d/iTuoo) = 0.163, where T denotes the period of the flow oscillation. This 
result is in a very close agreement with the experimental result by Williamson 
(1989), where 5 = 0.164 (see figure 3). 
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FIGURE 1 . Computational domain and grid system. The C-type mesh configura- 
tion is used for flow over a cylinder. The branch cut is located along the centerline 
of the cylinder wake. 



tUoo/d 

Figure 2. Time history of the streamwise velocity at several points (x = 2 Ad) 
behind the cylinder: , y = —0.21c?; , y = -0.52 d; , y = -l.ld; 

) V = —2. Id. Here, the position x = y = 0 corresponds to the center of the 
cylinder. 


It is interesting to note that the periodic behavior behind the cylinder was not 
found when an explicit interpolation method using adjacent points was used as 
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Re = Uocd/v 

FIGURE 3. Variation of 5 with Re: , the least square curve fit from 

Williamson (1989), 5 = A/Re + B + CRe , where A = —3.3265, B = 0.1816, C — 
1.6 x 10~ 4 ; •, the present calculation. 

the branch-cut boundary condition. It turns out that the use of the explicit inter- 
polation method along the branch cut indirectly prescribes zero pressure- gradient 
across the branch cut. Apparently, the absence of oscillating pressure force across 
the branch cut precludes vortex generation behind the cylinder. 

3. Future plans 

Future research plan includes two phases. 

Phase 1 : Laminar test cases will be computed and compared with experimental 
results. We will also compare laminar flow solutions of the NACA 0012 airfoil with 
those of a flat-plate airfoil for non- zero angles of attack in order to investigate the 
effect of the finite thickness of the airfoil on laminar separation characteristics. 

Phase 2: We will first simulate turbulent flow over an airfoil with zero angle of attack 
using the dynamic subgrid-scale model and compare the solution with experimental 
results and with those by Reynolds averaged turbulence closure models. The angle 
of attack for maximum lift (or for maximum lift-drag ratio) will be considered 
subsequently. 
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